import numpy as np
import pandas as pd
import plotly.express as px
import math
import plotly.io as pio
import plotly.graph_objects as go

# pio.renderers.default='notebook'
pd.options.plotting.backend = "plotly"

Utilities

# utils
# def plot_state_vs_time(df):    
#     fig = df.plot()
#     fig.update_layout(height=700)
#     fig.show()
    
def plot_state_vs_time(df):
    fig = go.Figure()

    fig.add_trace(go.Scatter(
        x=df.index,
        y=df.iloc[:, 0],
        legendgroup="group",  # this can be any string, not just "group"
        legendgrouptitle_text=df.columns.get_level_values(0)[2],
        name=df.columns.get_level_values(1)[2],
        mode="lines",
        marker=dict(color="blue", size=10)
    ))

    fig.add_trace(go.Scatter(
        x=df.index,
        y=df.iloc[:, 1],
        legendgroup="group",
        name=df.columns.get_level_values(1)[1],
        mode="lines",
        line=dict(color="red")
    ))
    fig.add_trace(go.Scatter(
        x=df.index,
        y=df.iloc[:, 2],
        legendgroup="group2",  # this can be any string, not just "group"
        legendgrouptitle_text=df.columns.get_level_values(0)[1],
        name=df.columns.get_level_values(1)[2],
        mode="lines",
        marker=dict(color="green", size=10)
    ))

    fig.add_trace(go.Scatter(
        x=df.index,
        y=df.iloc[:, 3],
        legendgroup="group2",
        name=df.columns.get_level_values(1)[1],
        mode="lines",
        line=dict(color="purple")
    ))
    
    fig.update_layout(height=700)
    fig.show()

Forward Euler Numerical Integrator

$\underline{X(t_{k + 1})} = \underline{X(t_{k})} + dt\underline{\dot{X(t_{k})}}$

def euler_integrator(x_0, t_0, t_n, dt, dx):
    x = x_0
    t = t_0
    
    yield (x, t)
    
    while t <= t_n:
        # eval x_k+1 using x_k and t_k
        x = x + dt * dx(x, t)
        
        # t increament to t_k+1
        t = t + dt
        
        yield (x, t)

Runga Kutta 4 Numerical Integrator

https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods

$\underline{X(t_{k + 1})} = \underline{X(t_{k})} + dt(k_1 + 2k_2 + 2k_3 + k_4)\frac{1}{6}$

where:
$f(x,t) = \underline{\dot{X(t_{k})}}$
$k_1 = f(\underline{x_k}, t_k)$
$k_2 = f(\underline{x_k} + dt\frac{k_1}{2}, t_k + \frac{dt}{2})$
$k_3 = f(\underline{x_k} + dt\frac{k_2}{2}, t_k + \frac{dt}{2})$
$k_4 = f(\underline{x_k} + dt k_3, t_k + dt)$

def rk4_integrator(x_0, t_0, t_n, dt, dx):
    x = x_0
    t = t_0
    
    yield (x, t)
    
    while t <= t_n:
        # eval x_k+1 using x_k and t_k
        k1 = dx(x, t)
        k2 = dx(x + dt * k1 / 2, t + dt / 2)
        k3 = dx(x + dt * k2 / 2, t + dt / 2)
        k4 = dx(x + dt * k3, t + dt)
        
        x = x + dt * (k1 + 2*k2 + 2*k3 + k4) * (1/6)
        
        # t increament to t_k+1
        t = t + dt
        
        yield (x, t)

Assignment Questions

Q2. b)

$\ddot{x}(t) = -2x(t) - 10\dot{x}(t)$

ans:

Let

$x(t) = x_1(t)$

$\dot{x}(t) = x_2(t)$

We define the state vector $\textbf{X} = [x_1, x_2]^T$

We can rewrite the equations given in the question as the following system of 2 first order O.D.E.s.

$\dot{x_1}(t) = x_2(t)$

$\dot{x_2}(t) = -2x_1(t) - 10x_2(t)$

A = np.array(
    [
        [0, 1],
        [-2, -10]
    ]
)

def dx(X, t = None):
    return A@X

X_0 = np.array(
    [
        [1],
        [0]
    ]
)

T = 20

# Euler
dt = 1e-3 # delta time (seconds)
state_time = []
state_integrator = euler_integrator(X_0, 0, T, dt, dx)
for x, t in state_integrator:
    state_time.append({("Euler", "X_1"): x[0][0], ("Euler", "X_2"): x[1][0], "Time": t})
e_state_time_df = pd.DataFrame(state_time).set_index("Time")

# RK4
dt = 1e-3 # delta time (seconds)
state_time = []
state_integrator = rk4_integrator(X_0, 0, T, dt, dx)
for x, t in state_integrator:
    state_time.append({("RK4", "X_1"): x[0][0], ("RK4", "X_2"): x[1][0], "Time": t})
r_state_time_df = pd.DataFrame(state_time).set_index("Time")

state_time_df = pd.concat([r_state_time_df, e_state_time_df])

plot_state_vs_time(state_time_df)

Q2. d)

$2\ddot{x}(t) - 5\dot{x}(t) + x(t) = 0$

ans:

Let

$x(t) = x_1(t)$

$\dot{x}(t) = x_2(t)$

We define the state vector $\textbf{X} = [x_1, x_2]^T$

We can rewrite the equations given in the question as the following system of 2 first order O.D.E.s.

$\dot{x_1}(t) = x_2(t)$

$\dot{x_2}(t) = \frac{5x_2(t) - x_1(t)}{2}$

A = np.array(
    [
        [0, 1],
        [-1/2, 5/2]
    ]
)

def dx(X, t = None):
    return A@X

X_0 = np.array(
    [
        [1],
        [0]
    ]
)

T = 20

# Euler
dt = 1e-3 # delta time (seconds)
state_time = []
state_integrator = euler_integrator(X_0, 0, T, dt, dx)
for x, t in state_integrator:
    state_time.append({("Euler", "X_1"): x[0][0], ("Euler", "X_2"): x[1][0], "Time": t})
e_state_time_df = pd.DataFrame(state_time).set_index("Time")

# RK4
dt = 1e-3 # delta time (seconds)
state_time = []
state_integrator = rk4_integrator(X_0, 0, T, dt, dx)
for x, t in state_integrator:
    state_time.append({("RK4", "X_1"): x[0][0], ("RK4", "X_2"): x[1][0], "Time": t})
r_state_time_df = pd.DataFrame(state_time).set_index("Time")

state_time_df = pd.concat([r_state_time_df, e_state_time_df])

plot_state_vs_time(state_time_df)

Q4.

$\ddot{x}(t) = 2y(t) - sin(4t^2x(t)), x(0) = 1, \dot{x}(0) = 2$

$\ddot{y}(t) = -2x(t) - \frac{1}{2t^2x(t)^2 + 3}, y(0)=1, \dot{y}(0) = 1$

ans:

Let

$x(t) = x_1(t)$

$\dot{x}(t) = x_2(t)$

$y(t) = y_1(t)$

$\dot{y}(t) = y_2(t)$

We define the state vector $\textbf{X} = [x_1, x_2, y_1, y_2]^T$

We can rewrite the equations given in the question as the following:

$\dot{\textbf{X}}(t) = f(\textbf{X}) = \begin{bmatrix}x_2 \\ 2y_1 - sin(4t^2x_1) \\ y_2 \\ -2x_1 - \frac{1}{2t^2x_2^2 + 3 } \end{bmatrix}$

def dx(X, t):
    x_1 = X[0][0]
    x_2 = X[1][0]
    y_1 = X[2][0]
    y_2 = X[3][0]
    
    return np.array([
        [x_2],
        [2 * y_1 - np.sin(4 * t**2 * x_1)],
        [y_2],
        [-2 * x_1 - 1 / (2 * t**2 * x_2 ** 2 + 3)]
    ])

X_0 = np.array(
    [
        [1],
        [2],
        [1],
        [0]
    ]
)

T = 2

# Euler
dt = 1e-3 # delta time (seconds)
state_time = []
state_integrator = euler_integrator(X_0, 0, T, dt, dx)
for x, t in state_integrator:
    state_time.append({("Euler", "X_1"): x[0][0], ("Euler", "X_2"): x[1][0], "Time": t})
e_state_time_df = pd.DataFrame(state_time).set_index("Time")

# RK4
dt = 1e-3 # delta time (seconds)
state_time = []
state_integrator = rk4_integrator(X_0, 0, T, dt, dx)
for x, t in state_integrator:
    state_time.append({("RK4", "X_1"): x[0][0], ("RK4", "X_2"): x[1][0], "Time": t})
r_state_time_df = pd.DataFrame(state_time).set_index("Time")

state_time_df = pd.concat([r_state_time_df, e_state_time_df])

plot_state_vs_time(state_time_df)